#include "Cell.h"

extern const vec<double> e[Q_DIM];
extern const double w[Q_DIM];

const double Cell::M[Q_DIM][Q_DIM] =	{{ 1, 1, 1, 1, 1, 1, 1, 1, 1},
										{-4,-1,-1,-1,-1, 2, 2, 2, 2},
										{ 4,-2,-2,-2,-2, 1, 1, 1, 1},
										{ 0, 1, 0,-1, 0, 1,-1,-1, 1},
										{ 0,-2, 0, 2, 0, 1,-1,-1, 1},
										{ 0, 0, 1, 0,-1, 1, 1,-1,-1},
										{ 0, 0,-2, 0, 2, 1, 1,-1,-1},
										{ 0, 1,-1, 1,-1, 0, 0, 0, 0},
										{ 0, 0, 0, 0, 0, 1,-1, 1,-1}}; //M[i][j]

const double Cell::Minv[Q_DIM][Q_DIM] ={{ 1./9., -1./9.,  1./9.,     0,      0,     0,      0,     0,    0},
										{ 1./9.,-1./36.,-1./18., 1./6., -1./6.,     0,      0, 1./4.,    0},
										{ 1./9.,-1./36.,-1./18.,     0,      0, 1./6., -1./6.,-1./4.,    0},
										{ 1./9.,-1./36.,-1./18.,-1./6.,  1./6.,     0,      0, 1./4.,    0},
										{ 1./9.,-1./36.,-1./18.,     0,      0,-1./6.,  1./6.,-1./4.,    0},
										{ 1./9., 1./18., 1./36., 1./6., 1./12., 1./6., 1./12.,     0, 1./4.},
										{ 1./9., 1./18., 1./36.,-1./6.,-1./12., 1./6., 1./12.,     0,-1./4.},
										{ 1./9., 1./18., 1./36.,-1./6.,-1./12.,-1./6.,-1./12.,     0, 1./4.},
										{ 1./9., 1./18., 1./36., 1./6., 1./12.,-1./6.,-1./12.,     0,-1./4.}}; //M^{-1}[i][j]

double Cell::meq[Q_DIM];
double Cell::m[Q_DIM];
double Cell::s[Q_DIM];

Cell::Cell(void) {
	isObstacle = false;
	for (int l=0; l<Q_DIM; l++) {
		f[l] = 0.0;
		m[l] = 0.0;
	}
}

void Cell::setS(double se, double seps, double schi, double sq, double snu) {
	s[iRHO] = 0.0;
	s[iE] = se;
	s[iEPSILON] = seps;
	s[iJX] = schi;
	s[iQX] = sq;
	s[iJY] = schi;
	s[iQY] = sq;
	s[iPXX] = snu;
	s[iPXY] = snu;
}

void Cell::updateRhoAndU(void) {
	double _rho = 0.0;
	double _jx = 0.0;
	double _jy = 0.0;
	for(int j=0; j<Q_DIM; j++) {
		_rho += M[iRHO][j]*f[j];
		_jx += M[iJX][j]*f[j];
		_jy += M[iJY][j]*f[j];
	}
	rho = _rho;
	if (rho > 0)
		u = {_jx/rho, _jy/rho};
	else
		u = {0.0, 0.0};
}

void Cell::calcMomenta() {
	for(int i=0; i<Q_DIM; i++) {
		m[i] = 0.0;
		for(int j=0; j<Q_DIM; j++)
			m[i] += M[i][j]*f[j];
	}
}

void Cell::calcDF() {
	for(int i=0; i<Q_DIM; i++) {
		f[i] = 0.0;
		for(int j=0; j<Q_DIM; j++)
			f[i] += Minv[i][j]*m[j];
	}
}

void Cell::calcEQ() {
	meq[iE] = 0.25*alpha2*m[iRHO] + (1./6.)*gamma2*(m[iJX]*m[iJX] + m[iJY]*m[iJY]);
	meq[iEPSILON] = 0.25*alpha3*m[iRHO] + (1./6.)*gamma4*(m[iJX]*m[iJX] + m[iJY]*m[iJY]);
	meq[iQX] = 0.5*c1*m[iJX];
	meq[iQY] = 0.5*c1*m[iJY];
	meq[iPXX] = 1.5*gamma1*(m[iJX]*m[iJX] - m[iJY]*m[iJY]);
	meq[iPXY] = 1.5*gamma3*m[iJX]*m[iJY];
}

void Cell::collide() {
	calcMomenta();
	calcEQ();
	for (int j=1; j<Q_DIM; j++)
		m[j] = m[j] - s[j]*(m[j] - meq[j]);
	calcDF();
	updateRhoAndU();
}

void Cell::collide(vec<double> _u) {
	calcMomenta();

	//s[iJX] = 1.0;
	//s[iJY] = 1.0;

	meq[iJX] = _u.x;
	meq[iJY] = _u.y;

	m[iJX] = m[iJX] - (1.0/0.908248290464)*(m[iJX] - meq[iJX]);
	m[iJY] = m[iJY] - (1.0/0.908248290464)*(m[iJY] - meq[iJY]);

	calcEQ();

	for (int j=1; j<Q_DIM; j++)
		if ((j != iJX) && (j != iJY))
			m[j] = m[j] - s[j]*(m[j] - meq[j]);

	calcDF();

	//s[iJX] = 0.0;
	//s[iJY] = 0.0;
}

void Cell::setU(vec<double>& _u, double dxux, double dyux, double dyuy, double dxuy, double _rho) {
	//equilibrio comun
	m[iRHO] = _rho;
	m[iJX] = _u.x*_rho;
	m[iJY] = _u.y*_rho;
	calcEQ();
	m[iE] = meq[iE];
	m[iEPSILON] = meq[iEPSILON];
	m[iQX] = meq[iQX];
	m[iQY] = meq[iQY];
	m[iPXX] = meq[iPXX];
	m[iPXY] = meq[iPXY];
	calcDF();
	updateRhoAndU();

	//término de no-equilibrio de Skordos según Lallemand y Luo
	/*double ux = _u.x;
	double uy = _u.y;
	m[iJX] = ux;
	m[iJY] = uy;
	m[iRHO] = 1.0;
	m[iE] = -2.0 + 3.0*(ux*ux + uy*uy);
	m[iEPSILON] = 1.0 - 3.0*(ux*ux + uy*uy);
	m[iQX] = -ux;
	m[iQY] = -uy;
	m[iPXX] = (ux*ux - uy*uy) - 2.0*(dxux - dyuy)/(3.0*s[iPXX]);
	m[iPXY] = ux*uy - (dyux + dxuy)/(3.0*s[iPXX]);
	calcDF();
	updateRhoAndU();*/

	//término de no-equilibrio de Skordos según Mei et al.
	/*m[iRHO] = _rho;
	m[iJX] = _u.x*_rho;
	m[iJY] = _u.y*_rho;
	calcEQ();
	m[iE] = meq[iE];
	m[iEPSILON] = meq[iEPSILON];
	m[iQX] = meq[iQX];
	m[iQY] = meq[iQY];
	m[iPXX] = meq[iPXX];
	m[iPXY] = meq[iPXY];
	m[iPXX] = m[iPXX] - 2.0*(dxux - dyuy)/(3.0*s[iPXX]);
	m[iPXY] = m[iPXY] - (dyux + dxuy)/(3.0*s[iPXY]);
	calcDF();
	updateRhoAndU();*/

	//término de no-equilibrio de Skordos según foro de palabos
	/*for (int l=0; l<Q_DIM; l++) {
		double _ex = e[l].x;
		double _ey = e[l].y;
		f[l] = f[l] - _rho*(1.0/s[iPXX])*3.0*w[l]*( pow(_ex,2)*dxux + pow(_ey,2)*dyuy + _ex*_ey*(dxuy + dyux) );
	}
	updateRhoAndU();*/
}

void Cell::addU(vec<double>& _u) {
	calcMomenta();
	m[iJX] += _u.x*m[iRHO];
	m[iJY] += _u.y*m[iRHO];
	calcDF();
}

void Cell::multU(double factor) {
	calcMomenta();
	m[iJX] *= factor;
	m[iJY] *= factor;
	calcDF();
}
